Forecasting confined spatiotemporal chaos with genetic algorithms 
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A technique to forecast spatiotemporal time series is presented, 
it uses a Proper Ortogonal or Karhunen-Loeve Decomposition to en- 
code large spatiotemporal data sets in a few time-series, and Genetic 
Algorithms to efficiently extract dynamical rules from the data. The 
method works very well for confined systems displaying spatiotem- 
poral chaos, as exemplified here by forecasting the evolution of the 
onedimensional complex Ginzburg-Landau equation in a finite do- 
main. 
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Nonlinear time-series analysis provides tools to identify dy- 
namical systems from measured data [[!]]. The approach has 
been greatly developed in the last years as a powerful alter- 
native to linear stochastic methods in the modeling of irreg- 
ular time-series and provides, under the assumption of de- 
terministic behavior, useful recipes for system control, noise 
reduction, and forecasting. Applications of these techniques 
to situations of spatiotemporal chaos, however, is still in its 
beginnings There are two main reasons for this: a) 

the large attractor dimensions of spatiotemporally chaotic sys- 
tems, increasing with system size, poses serious difficulties to 
the standard methods of delay embedding and attractor recon- 
struction; b) the right choice of variables is far from obvious: 
whereas the time evolution of an observable at a particular 
space point could be enough in some particular situations, de- 
caying space correlations, and propagation phenomena turn 
this to be a poorly performing choice in most cases. 

A very efficient method for time-series prediction using 
Genetic Algorithms (GA) has been recently proposed in [[|] 
for nonextended systems. Comparatively small data sets are 
enough to use this technique, which makes it competitive in 
facing difficulty a) i.e., prediction in the presence of attrac- 
tors of large dimension. In some cases, even non-trivial func- 
tional forms of dynamical systems generating the data can be 
unveiled [|[|. In this Letter we extend the GA approach to 
the forecasting of confined spatiotemporal chaos. By this we 
mean the situation in which chaotic dynamics in an extended 
system is strongly affected by the presence of boundaries. Our 
interest in this situation, somehow intermediate between low- 
dimensional chaos and homogeneous extensive chaos, arises 
from its relevance to real experimental situations and 
from recent work leading to theoretical understanding: 
the boundaries break translational symmetry and the result- 
ing phase rigidity restricts the shape of the chaotic fluctua- 
tions allowed. This manifests for example in the appearance 
of nontrivial average patterns [E|pl] and in inhomogeneities 



in other statistical characteristics Under these circum- 

stances the Empirical Orthogonal Functions (EOFs) JTo| , p~T] ] 
obtained from a Proper Ortogonal Decomposition (POD, also 
known as Karhunen-Loeve decomposition) provide an excel- 
lent basis for describing the system dynamics. They are differ- 
ent from simple Fourier modes and contain information (op- 
timal in a precise sense) on the broken traslational symmetry. 
The amplitudes of the most important EOFs will be the vari- 
ables chosen in response to difficulty b). 

We now describe more in detail our method for spatiotem- 
poral forecasting, in which the POD is used to encode the 
large spatiotemporal data set in a few time-series, and the GA 
approach is used to obtain the corresponding forecasts. Given 
a time series of spatial patterns U(x, n), where n = 1, N 
labels the temporal sequence and x the M spatial points in 
a d-dimensional mesh, the POD decomposes the fluctuations 
around the temporal mean u(x, n) = U (x, n) — (U (x, n)) n 
into modes ranked by their temporal variance. As a result, a 
set of spatial EOFs and associated temporal amplitude func- 
tions are obtained. The EOFs 0;(x) (i = 1, ...,M) are the 
(orthogonal) eigenfunctions of the covariance matrix of the 
data C(x, x') = (u(x, n)u(x', n)) n and are the spatial struc- 
tures statistically more representative of the fluctuations in the 
data set. Temporal amplitude functions ai(n), describing the 
dynamics of the system, are obtained from the modal decom- 
position it(x, n) = Xh=i a i{ n )4>i(' x -)- ^ on ly K < M of the 
EOFs (the ones containing the highest temporal variance as 
measured by the corresponding eigenvalues) are used in the 
reconstruction process, the set of reconstructed patterns 

K 

u K (x,n) = y"]at(ra)(fc(x) (1) 

i=l 

is still the best approximation one can obtain by linearly com- 
bining K arbitrary spatial patterns multiplied by K arbitrary 
amplitude functions [|10|], Even more, it has been shown for 
several chaotic and even turbulent confined systems JIo| , pl] ] 
that taking a few dominating modes K « M provides a 
good approximation to the complete data set. 

Forecasting of the amplitude functions is performed with a 
Genetic Algorithm. In general, GA's are computational meth- 
ods to solve optimization problems in which the optimal so- 
lution is searched iteratively with steps inspired in the Dar- 
winian processes of natural selection and survival of the fittest 
Jl2|]. Here the optimization problem to be solved is finding 
the empirical model best describing the data, that is, find- 
ing the optimum function Fi that minimizes the difference 
Ef = Xm=i { a i( n ) ~ ai(n)) 2 between the values ai(n) of 
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each time series and the corresponding estimator given by 

ai(n)=F i [ai(n-l),a i (n — 2),...,a i (n — D)] , (2) 

with D + 1 < n < N. Finding Fi,i = 1,...,M amounts 
to identify the dynamical system behind the data set. Once 
found, Eq. (Q) can be used to predict the future evolution of 
the system. If D is large enough, the existence of the exact 
Fi's is guaranteed by Takens theorem and its extensions [[!]], 
but a smaller D can give approximate dynamics Fi with al- 
ready a reasonably low error Ei. In addition, we are not look- 
ing for all the M estimators but only for the K associated to 
the dominant EOFs. In our approach, the time-series associ- 
ated to each EOF are modeled independently. More general 
multivariate estimators, with each Si possibly dependent on 
different ctj's, may in principle be used, but we restrict to the 
choice (||) for algorithmic simplicity. 

The power of the GA resides in that a huge functional space 
is explored in order to find an optimal Fi. Each possible Fi 
is a formula consisting in a combination of numerical con- 
stants, variables, and arithmetic operators. This combination 
is stored in the computer as a symbolic string. The only lim- 
itation to the allowed functional forms (besides the limitation 
to arithmetic operations) is the maximum allowed length of 
the symbolic string. The search procedure begins by ran- 
domly generating an initial population of potential estimators 
Fi that will be subjected to the evolutionary process. The evo- 
lution is carried out by selecting from the initial population 
the strongest individuals, i.e. the functions that best fit the 
data, giving a smaller Ei. In practice, only a temporal part of 
the data set is used in this step (the training set), whereas the 
rest of the data are used later for validating the efficiency of 
the prediction method (validating set). The strongest strings 
choose a mate for reproduction while the weaker strings dis- 
appear. 'Reproduction' consists in interchanging parts of the 
symbolic strings (the 'genetic material') between the two mat- 
ing individuals. As a result, a new generation of individu- 
als (which includes the original 'parent' string) is generated. 
The new population is then subjected to mutation processes 
that change, with low probability, small parts of the symbolic 
strings. The evolutionary steps are repeated with the new gen- 
eration, and the process is iterated until an optimum individual 
is finally found or after a fixed number of generations. Further 
details about the implementation of the algorithm can be con- 
sulted in [£|]. 

The formulae Fi are only optimized for predicting the value 
of cii (n) in terms of the D amplitudes immediately before in 
time. We call this 'one-step-ahead forecast'. One can in prin- 
ciple iterate the formulae to obtain successively predictions 
for cii(n +1), a.i(n + 2), etc. But this will normally lead 
to results rapidly diverging with respect to the correct values 
because of error accumulation and amplification [Q]. 

However, GA's can be designed specifically to forecast val- 
ues of the time series not necessarily in the immediate future. 
For example, finding the function Fj minimizing the error 
between the actual series and the estimator 



ai T {n) = Fj [a,i(n - T), a;(n - T - 1), ...,a,i(n - D)] 



(3) 



with D + 1 < n < N, allows direct prediction of a,i(N + T), 
that is prediction T-steps ahead, without iteration. 
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FIG. 1. Spatiotemporal evolutions of U(x,n), as given by the 
CGLE for q = 0.12 (left) and q = 0.16 (right). Black corresponds 
to U = 0, and lighter gray to high values of U. 

Numerical results. To illustrate the forecasting method 
we generate a data set from numerical simulation of a 
well-studied model equation displaying spatiotemporal chaos, 
the one-dimensional Complex Ginzburg-Landau equation 
(CGLE) , supplemented with Dirichlet boundary conditions 
at the ends of a finite interval [|lljj. It is convenient for our 
purposes to write it as 



d t A(x, t) = q 2 (l + a)d 2 x A + A - (1 + ip)A\A\< 



(4) 



where q, a, and (3 are real and positive and A(x, t) is a 
complex-valued field. We solve it in the interval [0, it] so that 
the boundary conditions read A(Q) = A(ir) = 0. By simple 
scaling of the spatial coordinate one sees that this is equiv- 
alent to rewriting the equation with q = 1, but solving it in 
a domain of size L = n/q. Thus the parameter q is equiva- 
lent to an inverse system size, and decreasing it is equivalent 
to increasing system size. Following [|TT]] we fix a ~ 4 and 
j3 = —4 [p"4|]. For q < 0.2 the system displays spatiotem- 
poral chaos for most of the initial conditions. Decreasing q 
one encounters the regime of confined spatiotemporal chaos 
we are interested in before approaching homogeneous exten- 
sive chaos at large system sizes (q — > 0) [Q. According to 
JTT|], the correlation dimension of the dynamical attractor for 
q = 0.14 is 9.08. We sample our simulation every r = 0.1 
time units and at spatial locations separated A = 7r/100 space 
units, and follow it for 80 time units (800 samples) after dis- 
carding the initial transient starting from random initial con- 
ditions (this sampling leads to N — 800 and M — 100). This 
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will be our 'training set' to be feed into the GA. The simula- 
tion is then continued for a few more time units, to provide 
the 'validation set' which is hidden to the GA. It is used later 
to check the accuracy of the predictions. 
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FIG. 2. The first two EOFs (a and b), and the corresponding time 
amplitude functions (c and d) from the training set at q = 0.16. 

We choose as the basic field to be forecasted the modulus 
U(x,n) — \A(x,t = nr)\ of the complex field. The algo- 
rithm seems to perform slightly better in forecasting the real 
or the imaginary parts of A, but we use U to show that the 
algorithm works well with nonlinear combinations of the ba- 
sic dynamical quantities. In Fig. ([!]) we show parts of typical 
spatiotemporal evolutions for q = 0.12 and q = 0.16. Clearly, 
reducing q decreases the spatial scales, as corresponding to an 
effectively larger system size, but also the complexity of the 
evolution is increased. In both cases it is clear that the motion 
of the dynamical structures is constrained by the presence of 
the walls, as corresponding to confined spatiotemporal chaos. 

We solve Eq. (|) for q = 0.18, 0.16, 0.14, 0.12 and perform 
the POD on the fluctuations u(x,n) of the modulus around its 
temporal mean value in the resulting data sets. The number of 
relevant EOFs (which we define to be those accounting for at 
least 99% of the data variance [ 1 1 1) are respectively 9, 11, 13, 
and 15. We note that this confirms the expected approximate 
linear scaling of the number of EOFs with increasing system 
size L (oc q^ 1 ) Jig], It is somehow surprising that this exten- 
sive scaling appears even when chaos is not homogeneous, but 
still influenced by the boundaries. This fact has been observed 
in other systems before For illustrative purposes, we 

show in Fig. ^ the two most relevant EOFs from our training 
set at q = 0.16, and the corresponding temporal amplitude 
functions. The chaotic character of these series is evident. 

We next apply the GA to each of the amplitude functions 
of the relevant EOFs. We use the following parameters for 
all the values of q: number of generations in the evolution- 
ary process 2000, number of individuals in each generation 
120, maximum number of symbols allowed for each symbolic 
string 20, number of delays in (^) or (||) D = 18. Tuning of 
these parameters for each particular value of q would improve 
forecasting, but would make comparisons more difficult. Pre- 
dictions for the field u(x, t) are then build up by reconstruc- 



tion according to ([[]) with K the number of relevant EOFs 
defined above. In Fig. (^J) we show the one-step-ahead fore- 
casted fields, more concretely the prediction for the first step 
beyond the training set, n = 801. It is compared with the ac- 
tual numerical pattern in the validation set, for q = 0.12, 0.14 
and q = 0.16, displaying an excellent performance. 
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FIG. 3. The forecasted moduli fields (dashed line) as compared to 
the real ones (solid) for one-step-ahead prediction for several values 
of q. 

We quantify the quality of the prediction in terms of the 
mean square error e q (n): 



4(n) = 



M 



1 

M ^ 



(u K (x=jA,n)-u{x = jA in )) , (5) 



where u K (x,n) is the predicted pattern reconstructed from 
Eq. (jl|) and u(x, n) is the actual pattern from the validation 
set. As stated before, GA's can be used to predict future values 
some time steps ahead, without the need of iterating the one- 
step-ahead predictor (which early becomes useless because of 
the expected exponential growth of errors). Figure (Q) shows 
e q (n) as a function of n for q = 0.16 calculated from: a) the 
one-step-ahead prediction formulae obtained from the training 
set, but applied to obtain the pattern at step n from the previ- 
ous D values in the validation set; b) iteration of the one-step- 
ahead formulae starting from the last D data in the training 
set; c) five-steps-ahead prediction from a formula of the type 
(j^) with T = 5, obtained by the GA in the training set, and 
used into the validation set. We see that the improvement in 
accuracy is notorious when iteration is avoided. We note that 
the errors in methods a) and c) remain bounded even when n is 
far from the values from which the prediction formulas were 
estimated (i.e. the training set n < 800). This confirms that 
the method is not simply fitting data, but rather it has really 
found approximate dynamical rules within the deterministic 
spatiotemporal series. 



3 




situations of confined spatiotemporal chaos as exemplified by 
the CGLE in a finite interval. We mention here that we are 
exploring the possibilities of the method for prediction from 
noisy natural data sets. Results obtained in forecasting Sea 
Surface Temperature patterns in an area of the Mediterranean 
Sea Jl7[ ] are encouraging. 

We acknowledge financial support from CICYT (MAR98- 
0840) and DGICYT (PB94-1 167). 



FIG. 4. Errors as a function of n in the validation set, for 
q = 0.16. Circles: one-step-ahead prediction. Diamonds: iter- 
ation of the one-step-ahead formulae starting from the training set 
(n < 800). Squares: five-steps-ahead prediction. 
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FIG. 5. Mean error for one-step-ahead prediction in the validation 
set as a function of q. 

Figure (^) displays the average error < e q >, which is the 
temporal average of e q (n) with n in the validation range dis- 
played in Fig. (Q), as a function of q (for one-step-ahead pre- 
diction). Despite we are including more EOFs in the recon- 
struction for decreasing q, the prediction error shows a ten- 
dency to increase. This is a consequence of the increase in 
complexity (and in attractor dimension) of the dynamics by 
the effective increase in system size (« q^ 1 ). Since we keep 
the number of delays D fixed, the embedding of the data set 
becomes more incomplete at smaller q and the prediction de- 
teriorates. In addition, for smaller q the confined or boundary 
influenced character of the spatiotemporal chaos in the sys- 
tem is lost and a description in terms of local structures will 
be certainly more efficient [^] . 

In summary, we have presented a method to forecast the 
evolution of spatially extended systems based in the combi- 
nation of POD and GA's. The method performs very well in 



[1] H.D.I. Abarbanel, Analysis of Observed Chaotic Data 
(Springer- Verlag, Berlin, 1996); H. Kantz and T. Schreiber, 
Nonlinear Time Series Analysis (Cambridge University Press, 
Cambridge, 1997). 
[2] U. Parlitz and G. Mayer-Kress, Phys. Rev. E 51, R2709 (1995); 

U. Parlitz and C. Merkwirth, Phys. Rev. Lett. 84, 1890 (2000). 
[3] S. 0rstavik and J. Stark, Phys. Lett. A 247, 145 (1998). 
[4] G.G. Szpiro, Phys. Rev. E 55, 2557 (1997). 
[5] V.K. Yadavalli, R.K. Dahule, S.S. Tambe, and B.D. Kulkarni, 

CHAOS 9, 789 (1999). 
[6] B.J. Gluckman, C.B. Arnold, and LP. Gollub, Phys. Rev. E 51, 
1128 (1995); L. Ning, Y. Hu, R.E. Ecke, and G. Ahlers, Phys. 
Rev. Lett. 71, 2216 (1993). 
[7] S.M. Zoldi, J. Liu, K.M.S. Bajaj, H.S. Greenside, and 
G. Ahlers, Phys. Re v. E 58, R6903 (1998); also in 



chao-dyn/9808006 



[8] V.M. Egufluz, P. Alstr0m, E. Hernandez-Garcfa, and O. Piro, 



Phys. Rev. E 59, 2822 (1999); also in 3hao-dyn/9805003 



[9] M. Meixner, S.M. Zoldi, S. Bose, and E. Scholl, Phys. Rev. E 
61, 1382 (2000). 

[10] P. Holmes, J.L. Lumley, and G. Beikooz,Turbulence, Coher- 
ent Structures, Dynamical Systems and Symmetry (Cambridge 
University Press, Cambridge, 1996); L. Sirovich, Q. Appl. 
Math. 45, 561 (1987). 

[11] L. Sirovich, Physica D 37, 126 (1989); J.D. Rodriguez and L. 
Sirovich, Physica D 43, 77 (1990); L. Sirovich, J.D. Rodriguez, 
and B. Knight, Physica D 43, 63 (1990). 

[12] J. H. Holland, Adaptation in natural and artificial systems (Uni- 
versity of Michigan Press, Ann Arbor, 1992). 

[13] A. Alvarez, A. Orfila, and J. Tintore (preprint, 2000). 

[14] The parameters Co and p of the notation of Sirovich et al. [JHJ] 
are related to ours by Co = 1/a, and p = —1/(3. Their time and 
space variables need to be multiplied by factors p and p/co, 
respectively, to obtain our time and space variables. 

[15] Our parameter regime leads to the so-called 'defect chaos' in the 
large-size limit. See B.I. Shraiman, A. Pumir, W. van Saarloos, 
PC. Hohenberg, H. Chate, and M. Holen, Physica D 57, 241 
(1992). 

[16] S. Zol di and H. Greenside, Phy s. Rev. Lett. 78, 1687 (1997); 
also in 



ghao-dyn/9610007 



[17] A. Alvarez, C. Lopez, M. Riera, E. Hern andez-Garcfa, and J. 
Tintore (preprint |chao-dyn/ 9 9 1 1 1 2|, 1999). 



4 



